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^ ■ Abstract 

For a sensor network, a tractable spatially-dependent node deployment model is presented with 
' the property that the density is inversely proportional to the sink distance. A stochastic model 

is formulated to examine message advancements under greedy routing in such a sensor network. 
The aim of this work is to demonstrate that an inhomogeneous Poisson process can be used to 
. model a sensor network with spatially-dependent node density. Elliptic integrals and asymptotic 

' approximations are used to describe the random behaviour of hops. Types of dependence that affect 

hop advancements are examined. We observe that the dependence between successive jumps in a 
' multihop path is captured by including only the previous forwarding node location. We include 

p I , a simple uncoordinated sleep scheme, and observe that the complexity of the model is reduced 

when enough nodes are asleep. All expressions involving multidimensional integrals are derived and 
evaluated with quasi-Monte Carlo integration methods based on Halton sequences and recently- 
developed lattice rules. An importance sampling function is derived to speed up the quasi-Monte 
Carlo methods. The ensuing results agree extremely well with simulations. 



^ ! 1 Introduction 

, Advancements in networking technologies are leading to sensor networks being a feasible and com- 

00 ' mon technology. However, a significant issue for sensor networks is developing routing methods that 

I can handle the dynamic topologies of sensor networks. A common approach in sensor networks and 

ff^ ■ ad hoc networks in general is geometric or position-based routing ^Wi [21] ■ The operation of these 

I algorithms is based on the assumption that each node knows its geographical location in relation to 

the main communication station, known as a sink, and the location of neighbouring nodes within 
its transmission radius. Geometric routing is often considered attractive because of its localized 
nature and scalability. A natural geometric routing approach is for the source node to forward a 
data message to the node that is geographically the closest to the sink, and repeat this step until the 
; ^ ■ message finally reaches the target sink. This approach serves as a greedy routing method in itself 

I or forms the basis for more intelligent routing methods in wireless ad hoc networks [3l [6l [121 HZ] • 

Sensor nodes are often randomly scattered over the sensor field. To conserve power-consumption, 
a subset of sensor nodes may randomly fall into a low energy-consuming sleep state, in which they 
cannot relay messages. The inherent randomness in sensor networks motivates the need for a suitable 
stochastic model to determine the ability of a routing scheme to successfully deliver data. In recent 
years, stochastic models and methods are being increasingly employed in analyzing communication 
networks. In particular, there is a special issue i51 on stochastic geometry and related fields applied 
to communication networks, as well as a two- volume monograph by Baccelli and Blaszczyszyn [Il[2]. 

The majority of this work, however, is under the assumption that sensor nodes are deployed ac- 
cording to a homogeneous Poisson process. Although such a tractable model might not capture the 
underlying node density variation, it can serve as a first approximation for studying network charac- 
teristics. We wish to extend the standard model by examining inhomogeneous or nonhomogeneous 
node deployment such that the node density is spatially-dependent. 

There are various reasons why inhomogeneous deployment models may be necessary. The de- 
ployment of sensor nodes depends on the environment and application of the sensor network. Con- 
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sequently, sensor nodes may need to be deployed more densely in important sensing regions. The 
obstacles in the network surroundings may prevent nodes being positioned in certain regions. In- 
consistent system parameters (such as battery lifetime) and interference can reduce the effective 
node density in certain regions. Furthermore, the nature of the actual node deployment influences 
the node density. For example, an aerial dispersal of nodes could result in the nodes obeying some 
type of diffusion process. Also, there may be network protocols that require positioning the nodes 
in certain regions, which lead to performance advantages. 

The design and deployment of sensor networks must address the problem of message collisions. 
The data messages in sensor networks converge towards the sink. Nodes closer to the sink need act 
as relay nodes more than nodes away from the sink. One possible solution to this problem is to 
deploy more nodes in these heavy traffic regions. Consequently, the node density would decay at 
some rate that is dependent on the distance to the sink. 

It is in this last setting that we wish to examine greedy routing. We propose an approach using 
a tractable mathematical model similar to the one developed in previous work [141 115] . Ideally, we 
want to offer a computationally quick and reliable way to obtain probabilistic descriptions of multi- 
hop paths in a sensor network with a simple stochastic sleep scheme. Moreover, we want to extend 
the model, analysis, and calculations methods from the homogeneous case to the inhomogeneous 
case, and demonstrate that the techniques are still valid. 

To achieve these goals, we analyze greedy routing in randomly deployed networks under the 
multihop situation. We propose a tractable spatially-dependent node density function, and analyze 
the resulting stochastic characteristics. Furthermore, we examine the influence of a simple stochastic 
sleep scheme. More specifically, we examine the effects of having a certain proportion of nodes 
awake at any given time. We derive resulting probability integral expressions, and demonstrate 
their feasible evaluation via quasi-Monte Carlo integration methods based on Halton sequences and 
recently-developed lattice rules. 

The work presented here is focused on the stochastic behaviour of multihop paths, the calcu- 
lation methods, and mathematically representing the effects of a sleep scheme. Overall, we show 
the application of this mathematical formulation in describing the stochastic behaviour of message 
delivery in sensor networks with inhomogeneous node deployment. Additionally, we demonstrate 
an implementable calculation procedure based on quasi-Monte Carlo methods for evaluating mul- 
tidimensional integrals. 

2 Background work 

The results presented here follow on from initial homogeneous Poisson model development and 
analysis |14j . which was later extended by evaluating resulting integrals and analyzing a simple 
sleep scheme [15) . Consequently, the majority of the formulation and techniques used here have 
been applied in the constant density setting. 

More generally, Ishizuka and Aida 11 performed some preliminary analysis of node placement 
approaches via simulations to gauge the fault tolerance of networks against random node breakdown 
and battery exhaustion. In particular, they assumed that individual nodes are scattered uniformly 
around the sink, and that the density decreased as the distance from the sink increased. Ishizuka 
and Aida examined two node density models: a Gaussian function and a simple power law (inversely 
proportional to the sink distance). They concluded that the simple power law placement resulted 
in good fault tolerance. 

A concept closely related to connectivity is the sensing coverage of a sensor network, which is 
the ability of a sensor network to successfully sense or cover the entire sensor field. Solutions to 
coverage problems have been based on coverage processes such as the Boolean model; see Hall [5] 
and Stoyan, Kendall and Mecke [3S] for more details. In the sensor network setting, a more recent 
example is that of Manohar et al. |5T] who used coverage processes to examine the coverage of a 
sensor network with an exponentially decreasing node density. 

The aforementioned work involving inhomogeneous node deployment cases did not examine the 
stochastic behaviour of any particular routing method. Furthermore, the work did not cover the 
effects that a sleep scheme has on stochastic dependencies in the routing model. 
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3 Mathematical model 



We assume that nodes communicate data radially, and that a node's transmission radius is a con- 
stant that clearly cuts off at some distance, which implies that a node can relay data only to another 
node when it is within the forwarding node's transmission radius. For numerical calculations and 
simulations, we rescale all lengths with respect to the transmission radius, and hence the transmis- 
sion radius is set to one. However, we denote the transmission radius r in ensuing calculations and 
equations for clarity and future extensions. 

We assume that nodes are scattered according to a two-dimensional Poisson process over a 
finite region, and that at any given time a random number of nodes are in sleep mode while the 
remaining are in their awake mode. Furthermore, we assume the awake node density, that is the 
average number of awake nodes per unit area, is a spatially-dependent function; that is, the nodes are 
scattered according to an inhomogeneous Poisson process [53]. We assume the nodes are scattered 
uniformly around the sink, but decreases in some way as the distance to the sink increases. Hence, 
the awake node density is a radial function of the form 

X{u) = \q{u) ue(0, cx)), (1) 

where u is the distance to the sink, and q{u) is a weighting function. The function A(u) can be 
interpreted as the mean number of awake nodes per infinitesimal area element. We refer to the 
constant A as the initial node density. 

Ideally, the q{u) function should be amenable to analytical and asymptotic methods while re- 
flecting a realistic node distribution. We study a node density that is inversely proportional to the 
sink distance 

l{u) = (2) 
u 

which is the model proposed by Ishizuka and Aida However, we chose this function also because 
we believe that the positioning of nodes will better accommodate the convergence of messages near 
the sink. 

We introduce a Poisson process mean measure A, which for a bounded Borel set i? C with 
area A is the density function integrated over the region B, thus in polar coordinates 

K{B) = ^ jj q{u)udude, (3) 

= \Q{B), (4) 

where Q is referred to as the reseated mean measure, and is used such that the notation is analogous 
to the results based on the constant density case when Q reduces to the area A of the region [T31 [TS] . 
We emphasize that all the Q-type expressions we consider in this work are derived in a similar 
manner to the corresponding area expressions under the homogeneous model by integrating the 
density over specific regions [HJ US] . 

Under a sleep model, the initial node density parameter can be written as A = pa where p 
is the probability of a node being awake and a is the underlying (that is, sleep and awake) node 
density parameter. The awake nodes form a thinned Poisson process. As in the homogeneous case, 
both the number of points kept and the number of points removed form random variables that are 
independent of each other [25] . Hence, in our model the number of awake nodes is independent of 
the number of asleep nodes. 

The number of awake nodes Nb located within some region i? is a inhomogeneous Poisson 
random variable with a probability mass function 

P(7Vs=n) = ^^^^e-^«(^). (5) 

4 Single hop analysis 

We introduce the parameter 7 to represent the distance between a forwarding node and the sink 
(refer to Fig. [1]). The sink distance parameter is set to 7 = ^ when the forwarding node is the source 
node. Often we shall present results such that the initial node density parameter A is a product of 
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the source-sink distance £ and some positive constant. 

For a forwarding node with a sink distance 7, let X^iu) C be its partial feasible region as a 
function of u. This region is formed by the intersection of two circles of radii r and u centered at 
the source node and the sink respectively (as illustrated in Fig. [Ij. The area of I^{u) is denoted 
by A^{u), and given by the integral 

A^{u) = 2 / wdedw, (6) 



where the sink angle function ^/'^(u) is defined as 



ip^{u) = arccos ( j . (7) 



2u7 

Furthermore, the mean measure of the region T^{u) is given by the integral 



A(I^(m)) = 2 / / X{w)wd0dw, (8) 



7 — r ^0 



which for the functions ((T|) and ([2]) reduces to 

A(I^(u)) = 2A / / dOdw. (9) 

J-l-r Jo 

Henceforth, we write 

A^(m) = A(X^(u)), Q-y{u) = Q{I^{u)), 

to refer to the mean measure and the rescaled mean measure respectively, the motivation of which 
will become apparent by the analogous results that follow. 

Let the random variable U be the sink distance of the forwarding node after a single message 
hop. The standard nearest neighbour [T3J [3S] method leads to the probability distribution 

1 M > 7 (10) 

M < 7 — r. 

There is a jump discontinuity in the distribution at u = 7 owing to the positive probability that 
no nodes lie within the feasible region. On the support where the distribution (fTUI) is absolutely 
continuous the probability density is defined by 



f{u)^\Q'{u)e-^Q-^''\ j-r<u<j, (11) 



'7 

where the derivative of the rescaled mean measure 

q;(m) - 2^^(u). 

For a node of sink distance 7, it follows from integral © that the mean measure as a function 
of u is calculated by the integral 

A^{u) = 2A / ip^{w)dw. (12) 

J -y — r 
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The corresponding indefinite integral has the analytic solution 

(7 + r) 



iE sin 



iF sin 



w 
7 + r 

w 
7 + r 



(7 - r) 

il + r) 
(7-r) 



(7 - r) 
(2r), 



(13) 
(14) 
(15) 



where i = ^/—l, and F and E denote incomplete elliptic integrals of the first and second kind 
defined in Legendre form as 



F{4>;k) 



Ei(t>; k) 



(1 - fc2 cos^ 



cos^ ( 



^1/2 



d0. 



(16) 



(17) 



A number of standard computer packages have pre-written routines for numerically evaluating 
elliptic integrals, however, care must be taken as their definitions vary depending on the package. 
The results presented here are obtained via our purposely-written elliptic integral functions based 
on Carlson symmetric elliptic integrals |5j. For more details see the papers by Carlson 4, 5 , which 
give algorithms that are readily implementable and can handle both real and complex values under 
specified variable and parameter regimes. 



The solution to the original integral ([T^ is 

+ 2iXE ( sin"^ 



u 



2iXE ^sin"^ 
2iXF ( sin"^ 



2iXF sin 



7 + r 
7 — r 
7 + r 
u 

7 + r 
7 — r 



(7- 



(7 + r 



(7 + r 



(7 + r 



7 + r / ' (7 — 



(7-r) 
(7-r) 
(2r) 
(2r). 



(18) 



The final result gives a real value, although in practice an insignificant imaginary part may arise from 
rounding errors. We found that the solution is obtained quickly owing to the speedy convergence 
of the Carlson algorithms. 

Let C = 7 — t7 be the distance advanced towards the sink when the originating node is at a 
distance 7. We let F denote the distribution of C, which leads to 

g-AQ^(7-c) o<c<r 

1 Or (19) 

c< 0. 



F^ic) 



Given that C is non-negative, the m-th moment expression 



E(C"^ 



m 



c'"-ip(C > c)dc, 



^m-lg-AQ^(7-c)^ 



follows; with the first moment given by 

E(C) = r- / e-^'^^^'^-'^^dc, 
Jo 

and the second moment given by 



L{C^) = r^ - 2 



(20) 



(21) 
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Figure 2: Comparison of h.i{u) asymptotic expressions (A = 3£, r = 1 and I — 10) 



4.1 Asymptotic results 

We derive an asymptotic approximation to the rescaled mean measure Q-y{u), which gives a tractable 
and accurate expression. The angle function ipj expanded around u = j — r gives 



w boiu - 7 + r)i/2 _ ^ + ^)3/2 ^ ^^(-^ _ ^ _^ ^^5/2^ 



where the expansion terms 
60 = 

bi = 



2r 



-,1/2 



7(7 - 
2r 
7(7 - r) 
2r 



li/2 r^2 



3r7 — 37 



21 



1I/2 r 



3r4 + 25r^j^ - 



lOr^-f + 30j^r - 57 



41 



.7(7-'^) 

follow. The asymptotic result for the rescaled mean measure 

'bo, 



1607-^(7 — r)^r'^ 



5 7 



(22) 

(23) 
(24) 
(25) 

(26) 



follows. 

The results of the second-order approximation are accurate for a unit radius (see Fig. [5]). Adding 
the third term only improves the results slightly. However, it may be needed for larger transmission 
radius models (as Fig. [3] suggests) . Conversely, the two-term expansion gives accurate results (see 
Fig. HI) when substituted into the sink distribution ([TU| . In fact, it appears that the expansion 
consisting of elementary functions can be used to give accurate results, which are clearly faster to 
evaluate than those based on elliptic integrals. However, we continue to use the exact solution of 
the integral (jl2[) . and later compare it and its approximation. 

We present asymptotic moment results for our spatially-dependent node density model. 



Theorem 4.1. For the mean measure 
moment 



E(C) 



provided "f > r, under greedy routing the first hop 

r(5/3) 



and the second hop moment 



2r 



r(5/3) 



r(7/3) 



(Ago) 



2/3 



(A,o)^/^' 



(27) 



(28) 
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Figure 3: Comparison of Ai{u) asymptotic expressions (A = 3^, r = 5 and £ 
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Figure 4: Comparison of Fi{u) asymptotic expressions {X — 3£, r — 1 and £ 
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as the initial node density 

A — >■ oo, 



where r(-) is the gamma function, and 



90=3 



2r 



7(7-0 



1/2 



Proof. Consider integrals of the form 



/(A) = / t'^e-^'^^'Ut, 

J a 

where Q{t) > Q{a) for all t e {a,b), and asymptotically 

Q(t) ^ qo{t — a)'^ as t ^ a. 
The generalized Laplace's method [2Bj page 58] applied to this integral gives the asymptotic result 

T(t) 

tJ'iMay 

where 

2(fc + 1) 

For Laplace's method, Q{a) needs to be the minimum on the integral interval. Hence, the change 
of variable t = u — ^ + r = r — c with a slight abuse of notation leads to the area function expansion 

Q^{t)^qot^^^ + qit'^/^ + 0{t^/^) as i^O, (29) 

where 

46o 3 

The first moment result ([27]) follows, and the change of variable applied to the second moment 
equation gives 

E{C^) = - 2 Hr ~ t)e-^Q^^'Ut, 
Jq 

which leads to the second result ([25]) . □ 



4.2 Sink dependence 

The node intensity function is clearly dependent on the forwarding node sink distance. Compar- 
ing the hop distributions (in Fig. ^ reveals that a message is relayed farther in a single hop if 
the forwarding node is closer to the sink. Intuitively, hops increase stochastically as the message 
approaches the sink as more potential forwarding nodes are available in the forwarding region. 
Geometrically, we observe that the integral kernel in the mean measure equation (|12p is the angle 
function ipj, which decreases as 7 increases; that is 

-071 (") < "072 (^); 71 > 72 > r. (30) 

Conversely, tpj and, hence, the indefinite integral increases as the sink distance decreases. This 
dependence on the sink distance of the forwarding node is simply referred to as the sink dependence. 

The influence of the sink dependence can be observed by comparing the difference in two hop 
distributions with different sink distances 71 and 72. In previous work [14j . the hop distribution 
dependence on the sink distance under the homogeneous model was examined by a KuUback- 
Leibler analysis. Subsequently, the Kullback-Leibler divergence [TB] applied to the mixed discrete- 
continuous hop distribution gives 

^^(71,72) = / /72(c) log 
Jo 



fj2 (C) 
./71(C). 



dc-f F^3(0+)log 



^"72(0^ 

.^71 (0^ 



(31) 
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Figure 5: Numerical and asymptotic results for first moment (r = 1 and £ 
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Figure 7: Kullback-Leibler analysis of two hop distributions for (A = 3^, r = 1 and £ — 10). 



where the routing void (that is, no nodes in the feasible region) probability 

FJQ+) = limF^(e), 

The Kullback-Leibler divergence is non- negative and is zero for identical distributions jl6| . 

We calculated the integral in equation (pij) numerically to observe how the hop distribution 
is influenced when we set 71 = £ and vary 72 — 7. The comparison reveals that D(£,j) is high 
near the sink and decreases as 7 increases (refer to Fig [7]). This is the same expected behaviour 
as observed under the homogeneous model [H], however, D{£, 7) decays relatively slowly under the 
inhomogeneous model. This contrasts starkly with the homogeneous model where D{£, 7) is mostly 
zero away from the sink, and only increases significantly near the sink |14| . Furthermore, Z)(€, 7) 
differs significantly for two different £; the forwarding node sink distance cannot differ greatly when 
calculating hop distributions. This distribution behaviour eliminates the possibility of using renewal 
processes to model message advancement over multihop routes [37J [T3] . 

5 Multihop analysis 

For a multihop analysis, indexing for the random variables U and C is introduced. We set Uq = £ 
and Ui ~ U such that the i-th hop advancement is Ci — Ui^i — Ui. Each Ci depends on the 
forwarding node's sink distance Ui-i\ thus, the sink distance of the source node clearly affects the 
first hop and subsequent hops. As noted in the previous section, comparing the hop distributions 
demonstrates that each Ci is stochastically dominated by C^+i. That is, for i > 0, we have the 
stochastic ordering 

nC^+l > c) > P(a > c), c e (0, r). (32) 

This inequality is the opposite to the equivalent result under the homogeneous model as noted by 
Zorzi and Rao |17] and Keeler and Taylor |14) . 

5.1 Path dependence 

Let the random variable Qi be the angle between the i-th node and the previous node in relation 
to the sink. We assign the point Xi = {Ui, Qi) to the j-th forwarding node. The source (or zcroth) 
node corresponds to the point Xq = {£,0). A message travels i hops along a path that corresponds 
to a sequence of random points Xi = {Xq, Xi, . . . , Xi). 

Let Ii{ui+i) C be the feasible region of the i-th forwarding node as a function of Ui+i under 
the independent model. After the first hop, the nature of greedy routing implies that another 
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Figure 8: No awake nodes in the intersection region Ii \ Iq during the first message relay. 



dependence arises in the distribution of C/2, which was observed under the homogeneous model 
[l4l [15] . If a forwarding node is chosen, then there are no other nodes in the source feasible region 
closer to the sink. Hence, the intersection of the feasible regions of the source and the first node 
(that is, Xi nio in Fig- 13) has no awake nodes. This implies that U2 is dependent on both Ux 
and 61, the angle between the first node and the source node in relation to the sink. We call this 
dependence in both the sink distance and the sink angle the path dependence, and the hop model 
that includes both the path and the sink dependence simply the dependent model. Conversely, the 
independent model only includes the sink dependence. 

Under a sleep scheme, if nodes were asleep during the previous message relay, it is possible for 
recently awoken nodes to be present in this region during the current message relay. To perform 
the initial analysis we assume no sleep scheme exists (by setting p = 1). After our analysis, we 
include a simple sleep scheme and examine how varying p and a (while fixing A) affects the path 
dependence. 

The distribution of Ui+i under the independent model is dependent only on the sink distance 
of the current forwarding node, hence we write 

F, (u,+i) = P/(C/,+i < u,+l\U^ = u,), (33) 

while under the dependent model the distribution is dependent on the message path, and so we 
write 

G, (m,+i) = Pd(C/,+i < u,+i|X, = X,), (34) 

where the subscripts / and D indicate probability measures under the two models (the subscript is 
dropped if a result applies to both models). We denote the rescaled mean measure of the feasible 
region under the independent and dependent models respectively as Qi{ui+i) and Qi{uiJ^i). 

The rescaled mean measure under the independent model is always given by the original equation 
(jlSp . and hence, the distribution and probability density of J7i+i are obtained by setting 7 = in 
equations ((TU]) and (TTT]) . Under the dependent model, if the rescaled mean measure is given after i 
hops, the sink distribution 

{l_ f,-\Q,{n,+i) Ui-r<u<Ui 
1 u>u, (35) 

u < Ui ~ r, 

immediately follows, and where it is absolutely continuous its probability density 

g,{u,+i) = Ag:(?.,+i)e-^'5,K+i)^ (36) 

also follows. 

The set representing the feasible region under the dependent model follows by excluding the 
intersections of previous feasible regions, namely 

Viiui+i) = Ii{ui+i) \ U^~o"%-(uj+i). 

Under the dependent model, it is possible to calculate the rescaled mean over the feasible region 
after one hop; see appendix for details. If i > 2, we approximate the feasible region under the 
dependent model 

T)i{ui+i) « Ii(ui+i) \Ii_i(uj+i), (37) 
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where we refer to this approximation as the one-hop modeL Resuhs under the homogeneous model 
revealed that this approximation sufficiently captures the path dependence [15[ . Consequently, only 
the location of the previous node is needed to calculate the rescaled mean under the dependent 
model. 

The path dependence depends on the precise locations of the previous nodes while the sink 
dependence which only depends on their distance from the sink. To describe the random behaviour 
of message hops we need the joint density of Ui and Oj. Under our inhomogeneous Poisson model, 
the joint probability density is 

g^iu,+l,e,+l) - Ap.K+i,0,+i)e-^'5.K+i)^ (38) 

where the derivation is analogous to that of the homogeneous model [T5]; see appendix for details. 
The spatially-dependent initial density Ap. (wi+i, Oi+i) = Alp,. (ui^\, ^i+i) and the indicator function 
of the dependent feasible region 

Ip,K+i,0,;+i)-| otherwise, 

where angular coordinate Q^i is the angle between the source node and the «-th forwarding node in 
relation to the sink. 



6 Sleep model 

We outline a simple sleep scheme that has been analyzed under the homogeneous model [M]. We 
assume that the probability that a node is awake on each hop is p and the event that a node is awake 
during a transmission attempt is independent of the event that it is awake at other transmission 
attempts. Consequently, the intersection region {X\ nio in Fig. [5]) has a thinned initial node density 
(1 — p)A. The rest of the feasible region X\ \ Iq has an initial density A. It follows that the initial 
node density function after one message hop is given by 

Xv^{u2,02) = A [Ii,\i„(u2,02) + (1 -p)Iiinio("2,^^2)] , (39) 

where the superscripts denote the indicator functions of the disjoint regions. 

In the limit as p approaches zero and a approaches infinity with A = pa held constant, the 
locations of the nodes after each hop is re-sampled, thus, completely removing the path dependence 
and allowing each forwarding node to be treated like a source node. Its effects on the node density 
have been explored more thoroughly under the homogeneous model [15^. 

To calculate the Poisson mean measure over the region Vi, the above node density function is 
integrated over the domain [3S] leading to 

Ax>,(ui+i)= / / Xv,{wi+i,6i+i)q{wi+i)wi+id9i+idwi+i, (40) 

and we define 

Q-Di{ut+i) = . (41) 

The joint probability density of Ui and Oi, that is 

g^{u,+u9,+i) = Ap.(?/,+i,&,+i)e-^'5.,.K+i). (42) 



7 Multihop distribution 

7.1 Hop advancements 

We initially formulate this problem with the sink distance variabe since greedy routing is naturally 
based on it. However, hop advancement is a more intuitive variable in describing message progress 
over a multihop path. We adopt similar notation used for the sink distance random variables, hence 
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F and G denote the distributions under the two models. The complement of the sink distribution 
yields the hop distribution under both the independent and dependent models; the latter being 

( e-^'3'("'-^'+i) < C^+l < T 
G«(c,+i) = <^ 1 Q+1 > r (43) 

I Q+i < 0. 

and its probability density defined on the absolutely continuous part of the support 

5,(c,+i) = \Q\{u, - c,+i)e-^'3.K-c,+i)^ (44) 

where a simple sum relates the hop and sink distance variables 

i 

u,^l-Y,c,. (45) 



7.2 Distribution of Z 



n 



Let the random variable Z„ represent the distance advanced by a message in n hops 

n 

^„ = ^C,. (46) 

i=l 

To calculate the distribution of Z„ we use the joint probability density of the random variables Ci 
to Cn and 81 to 0„, that is 

n 

5(n-i)(ci, c„, 01, ., On) = W Ap._, (u,_i - c„ 0,)e-^'3.-iK-i-c,)^ (47) 

i=l 

which is defined for Ci G (0, r]. The derivation of the joint probability density is similar to that under 
the homogeneous mode [15j . which we have adapted and included in the appendix for completeness. 
It follows that the distribution of message advancement after n hops is expressed by 

(■min(z — ci ,r) ri>l{c2) 
AZn <Z)^ I dci I dOi / dC2 ... (48) 

'0+ 01(C2) 

(ci,.,c„,0i,.,0„)d0„ (49) 
+ ^d{Ci = 0) + Pz5((Zi < z) n (C2 = 0)) + . . . (50) 

+ Pi,((z„_i <z)n(c„ = o)). (51) 

where 

^i{ci+i) = ipui{ut - Ci+i), (52) 

denotes the maximum angle for a sink distance given by the sink angle function ([7])- The distribution 
of Zn under a sleep scheme is obtained by substituting the product of the joint probability densities 

(sa). 

The integral explicitly shown in the expression of V]j(Zn < z) is the distribution of Zn condi- 
tioned on the event that all hops Ci advance some positive distance. We will often refer to this 
integral simply as the conditional distribution of Zn, and denote it by 

P(Z„ < z\+) = P{Zn < z\Ci > 0, . . . , C„ > 0). (53) 

Under the independent model the joint probability density is not a function of any of the variables 
9i to 9n- Hence, the equivalent integral can be analytically integrated over the sink angle domains 



pmin(z. 




-0o(ci) 




do, J 




'0+ 




-■4>o{ci) 


nmin{z- 










dCn 


'0+ 







13 



[14[ , thus giving a simplified expression in the form of hop probabihty densities 

nmin{z^r) nmin{z — ci ,r) 

F/(Z„ < z) - / dci dc2... (54) 

/'inin(j2 — Ci.r) 

/o(ci) . . . /„_i(c„)rfc„ (55) 

0+ 

+ P/(Ci = 0) + P/((Zi < z) n (Ca = 0)) + . . . (56) 
+ P/((^„-i <z)n(C„ = 0)). (57) 

Under the dependent model, the probabihty of a message reaching a routing void after advancing 
i hops 

Voia+i - 0\X., = X,) = e-^Q.K)^ (58) 

follows. The routing void probability leads to the distribution of Z„ conditioned on the event that 
the message meets a routing void on the last hop 

min(z,r) /•4>o(ci) 



^D{{Zn <Z)C^ {Cn+1 = 0)) = / rfci / dOi 



0+ J'il>o{ci 
min(2 — Ci,r) 

dCn-- ■ 

+ 

5(„_i)(ci,.,c„,6'i,.,6'„)P_d(C„+i = 0)d6i„. (59) 

Consider the event when a message does not advance, hence Xi = X^+i. If there is a sleep 
scheme, a forwarding node can benefit by making multiple relay attempts. The number of different 
possible events results in the integral expressions needed to describe such a model soon growing 
to be intractable. We restrict the integrals by assuming that a message executes only one relay 
attempt. We point out that under the homogeneous model [15 , stochastic 'rules of thumb' have 
been proposed for how many re-attempts forwarding nodes should make before considering other 
options (such as message backtracking). Under the inhomogeneous model, the equivalent results 
can obtained by appropriately replacing the area terms with the corresponding integrals from the 
rescaled mean measures. 



7.3 Number of hops 

Let the random variable N represent the total number of hops required for a message to reach the 
sink. The distribution of N gives a different perspective of the performance of a routing method 
in a sensor network. The random variables N and Z„ are connected by a simple result [13], which 
has been used to calculate the distribution of N from the distribution of Z„ in the homogeneous 
setting [TS]. We use this relation to calculate the equivalent results under our spatially-dependent 
node density model. 



8 Integration methods 

To calculate the distributions of Z„ under the dependent model, a 2n-fold integral needs to 
be evaluated. For low n, traditional numerical integration schemes can be used. Unfortunately, 
integration by these methods is too slow at higher hop numbers, which motivates us to employ 
quasi-Monte Carlo methods. 

8.1 Quasi-Monte Carlo 

The integration description that follows is similar to the more detailed account [T31 Chap. 4] where 
quasi-Monte Carlo methods were used to evaluate similar integrals under the homogeneous model. 
In recent years, these integration methods have gained much interest owing to their speed and 
accuracy in evaluating high dimensional integrals. 
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Quasi-Monte Carlo methods are based on purely deterministic sequences of quasi-random num- 
bers such as Halton ^Oj and Sobol [22^ sequences. Mathematically, quasi-random sequences have 
low-discrepancy j20l 123) . Informally, such sequences appear 'less random' than sequences produced 
by regular pseudo-random number generators as they occur more evenly spaced apart. 

Previous numerical work iJJj has led us to use leaped Halton sequences to calculate integrals 
mostly. Often integrals were also calculated via regular Monte Carlo methods to check the quasi- 
random approach. We also give some results based on so-called lattice rules, which lead to specific 
cases of quasi-random sequences. These rules can produce well-behaving quasi-random sequences 
and have been a research focus in recent years owing to their ability to counter the curse dimen- 
sionality jl8j . The points arising from lattice rules are used in a similar manner to the quasi-Monte 
Carlo approach. There are many suggestions for lattice rules, but we lightly examine only one 
known as rank-1 lattice rule, which over the unit hyper-cube gives the integral estimate 



where the generating vector z G Z*, 71 is the number of function samples, and the braces give the 
fractional part in [0, 1). 

Given a generating vector, quasi-random sequences based on such a lattice rule can be clearly 
produced in an exceedingly fast manner. The way to quickly evaluate an integral is by choosing a 
suitable generating vector. However, the drawback is that lattice rules are based on input parameters 
known as weights, which depend on the nature of the function. In particular, the weights depend 
on how the function varies with respect to all its variables. Furthermore, the choice of some lattice 
rules require the total number of function samples before the integral calculation starts. This 
differs from regular quasi-Monte Carlo methods in which the integral estimate can be calculated 
continually until a sufhcient number of function samples has been taken. 

Under lattice rules, the number of function samples also influences the choice of the generating 
vector. Consequently, using the most suitable generating vector may not be a simple task as it 
involves analyzing the function of interest. However, a thorough analysis of which weights and 
quasi-random sequences are the most suitable in this setting is beyond the scope of this work. We 
simply give some complementary results based on the rank-1 lattice rule, and leave the analysis as 
a future task. 

The research field of lattice rules has a relatively short history, and new lattice rules are being 
developed continually. For more information, we refer the reader to the introductory piece by Kuo 
and Soan |18j , and an example of lattice rules used in a financial setting ^7 . Suitable lattice rules 
may offer a substantially faster way of evaluating the high dimensional integrals presented here and 
in previous work |15] . However, we focus on reducing error regardless of the chosen sequences by 
employing importance sampling, which has been done under the homogeneous model |15j . 

8.2 Importance sampling 

To reduce the variance of the integral estimate we employ importance sampling; that is, suitably 
generate Ci values to sample the function in key regions. Based on previous work [15j . we derive an 
importance sampling function which is similar in form to that of the homogeneous case. We recall 
the expansion of the Q function (|26p in which we use only the first term, hence 




(60) 



as 



7 + r)3/2 + 0(u)5/2 



(61) 



M — >■ 7 — r. 



where 



3 



A change of variable c = 7 — it leads to the function 



Q^(c) = go(r-c)3/2, 



which leads to an approximate solution for the hop distribution 




~ e 



AQ^(c) 
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The quantity ^^(O) is subtracted from the above expression and the result is divided by AF-y 



F^(cniax) — ^7(0)1 to obtain an importance samphng function 



1 



< C < Cn 



(62) 



where Cmax is the largest hop value. The subtracting of the routing void term has negligible effect 
for sufficiently large A. The corresponding derivative needs to integrate to one, hence the rescaling 
step. The derivative 

3A 



/7(c) 



2AF. 



1/2 -AQ^(c) 



(63) 



exists on the same domain as the sampling function. The inverse of the sampling function is given 

by 



In 



tAF^+F^{0) 



0<t<r. 



(64) 



Substituting a random variable from a uniform distribution, say T ~ C/(0,r), into the inverse of 
the importance sampling function gives a random variable C adhering to the importance sampling 
distribution (15^ . 

After each hop is generated, the functions ((63|) and (|64| should be updated by setting 7 to the 
current sink distance for each sample. This step was not necessary under the homogeneous model. 
However, under this inhomogeneous model each hop distribution varies more with respect to the 
sink distance. Furthermore, for high A, the importance sampling step should improve as it it based 
on the independent model, and at high node density the intersection regions grow stochastically 
smaller [HIIS]. 



9 Simulation 

We used routing simulations to verify our stochastic model and calculations. Given a source node 
sink distance £, message relaying is simulated in a circular sensor field Ci C of radius £ with 
the sink located at the origin. The fact that messages only advance towards the sink under greedy 
routings implies that sensor field edges do not infiuence the message routing. The total number of 
nodes per simulation is a Poisson random variable with the parameter 

K{Ce) = A / / q{u)udude 
Jo Jo 

= 2Xn£. 

In simulation, a node is assigned two independent random variables Qs and i?s corresponding to 
their polar coordinates in relation to the sink. To simulate node deployment such that the nodes 
adhere to the spatially-dependent node density ([2]) , both random variables are uniformly distributed 

P(es<0)=-, 0e[O,Tr], 

TT 

nRs<r)^j, re[0,^]. 

10 Results 

All the numerical integration and routing simulations were performed in Matlab on a standard ma- 
chine. The built-in Halton sequence generator was employed for the quasi-Monte Carlo integration. 
For a given value of z, it took between 10^ to 10^ points to obtain a quasi-Monte Carlo estimate 
of the conditional distribution P(Z„ < z\+). The importance sampling step improved the rate of 
convergence, particularly for high A. The actual number of points depends on the number of hops; 
more hops require more points. All calculations took no longer than an hour to complete, and 
usually considerably less. 
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Figure 9: Analytic and empirical distributions of the sink distance after one hop (p = 1,A = 3€, 
r = 1 and £ = 10). 



We compared the dependent model to routing simulations of various ensemble sizes with and 
without a blinking sleep scheme. Generally, between 10'^ to 10^ routing simulations were required. 
Similar to the integration process, higher hop numbers required more simulations to give converged 
results. The routing simulations are based on the same assumptions made in the mathematical 
model. 

We observed under the inhomogeneous model that more function samples and routing simu- 
lations are needed to give similarly converged results compared to those obtained under the ho- 
mogeneous model. An extensive empirical investigation is needed to see how fast quasi-Monte 
Carlo methods are compared to routing simulations. Also, more empirical and theoretical evidence 
is needed to elucidate the advantages and disadvantages of calculating probabilistic behaviour of 
greedy routing via our model. 

For a single hop, the mathematical model results are practically identical to the simulation 
results (see Fig. ^ . This implies that the model adequately captures the spatially-dependent node 
density function, and also vindicates the use of the elliptic integral functions. 

The lattice rule approach was only used to calculate one set of results (see Fig. [TT]). We used 
generator vectors based on fixed lattice rules for 2^° function sample points and equal weights [T7] . 
We applied ten random shifts to the lattice rules; see Giles et al. [7] for an example. This approach 
generally performed well, however, under importance sampling sometimes erratic results arose. This 
is possibly due to an over bias in function sampling or an unknown numerical artefact. A more 
thorough examination is needed, but we believe that the preliminary results using lattice rules are 
promising. 

Despite the accuracy of elliptic integral functions, we found that the results based on them could 
be replaced with results that used the three-term approximations and (pS)) with no discernible 
loss of accuracy (see the plots in Fig. 1111) . Evaluating the approximations is considerably faster as 
the expressions only involve elementary functions. Consequently, the remaining results are based 
on these approximations (Fig. [TUto Fig. [H]) . 

We calculated the distribution of N, and included a simple sleep scheme to see its influence on 
the path dependence. We compared the independent model to the dependent model by varying 
p (and accordingly a, hence holding A constant). The difference between the two models has an 
accumulative effect on N, hence it serves as a good indicator of path dependence. 

Under the independent model, the distribution P(A^ < n) gives greater values for each n than 
the equivalent result under the dependent model. Also, the path dependence clearly lessens as p 
approaches zero (compare Fig. [T^ and Fig. [T51) . These results are analogous to those under the 
homogeneous density model |15] . 

For large A, the difference between the two models is less. As is the case for the homogeneous. 



Analytic 

Simulation 
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Figure 10: Results of P(Z2 < z\+) via quasi-Monte Carlo integration and simulations {p = 1,A = 3£, 
r = 1 and ^ = 10). 




Figure 11: Results of P(Z2 ^ z\+) via lattice rules integration (based on analytic and asymptotic 
expressions) and simulations {p — l,X — 3£, r — 1 and £ = 10). 
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Figure 12: Independent and dependent model results oiF{N < n) compared to simulations (p = 1, 
A = 2£, r = 1 and £ = 10). 
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Figure 13: Independent and dependent model results of P(A^ < n) compared to simulations (p = 0.1, 
A = 2£, r = 1 and £ = 10). 
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Figure 14: Independent and dependent model results ol P(iV < n) compared to simulations {p — 1, 
A = 3^, r = 1 and £ = 10). 



a larger density results in the next forwarding node being closer to the sink, thus reducing the 
intersection of the feasible regions and lessening the path dependence. Moreover, under the inho- 
mogeneous model hops grows stochastically larger due to the increasing node density. Thus, we 
believe that the path dependence continues to decrease stochastically as the message approaches 
the sink. 

In conclusion, the results reveal that the dependent model clearly captures the path dependence. 
The resulting expressions, both involving the elliptic integrals and the asymptotic expansions, give 
results that closely agree with simulations. We believe more numerical investigation is needed to 
choose the most appropriate quasi-random sequences in this setting. 



11 Future work 

In realistic settings, the constant node density assumption may often not be appropriate. Under 
the spatially-dependent model, the density function was chosen such that it was simple enough 
for analytic and asymptotic means, while still being a plausible node placement scenario. Other 
suggestions exist such as the node density decaying exponentially or according to some power of 
the sink distance. 

Furthermore, under our model the sink was located at the maximum of the node density. Placing 
the sink at an arbitrary point in the sensor field results in the node density being dependent on the 
sink angle. This is an additional increase in the complexity of the density function. Moreover, the 
angle of an individual node would not be a uniformly distributed random variable, thus importance 
sampling might be needed when integrating over the angle domains. Consequently, these suggestions 
may result in analytic and asymptotic mean measures that can be used to model more realistic node 
deployment models. 

Further investigation of the asymptotic approximations are needed. The approximation may 
break down when the radius is large compared to the sink distance. For a constant radius model, 
the lengths can always be rescaled with respect to the transmission radius. However, this may not 
be possible for a randomly varying radius model. We stress that including random transmission 
radii into our model would be an interesting and realistic model extension in itself. 

We presented some integration results based on lattice rules. Despite these methods giving 
mostly agreeable results, further investigation is needed to gauge which lattice rules and quasi- 
random sequences are the most suitable for the integrals that arise from our model. This investi- 
gation would be both analytic and numerical in nature. This work may lead to our mathematical 
model considerably outperforming regular routing simulations. 

Finally, an attractive feature of regular Monte Carlo methods is that the error is obtained by 
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estimating the variance of the integral. Conversely, quasi-Monte Carlo methods lack a practical 
way of estimating the error, despite them generally have a faster convergence rate. The idea of 
'randomized' quasi-Monte Carlo methods seeks to combine the advantages of both approaches. 
Consequently, a future task lies in investigating these hybrid methods in evaluating hop integrals. 

12 Conclusion 

We presented a tractable inhomogeneous spatially-dependent density model. Inspired by previous 
work, we developed and examined a greedy routing model that incorporates both sink and path 
dependence. Moreover, the spatially-dependent density model verified that the formulation of the 
homogeneous model can be extended to an inhomogeneous case. This model is an alternative means 
of ascertaining the stochastic characteristics of greedy routing in sensor and ad hoc networks. 

We used asymptotic methods to derive accurate approximations to hop length moments and 
the mean measure for our spatially-dependent node density model. We uscid quasi-Monte Carlo 
methods and recently-developed lattice riiles coiipled with importance sampling to estimate the 
resulting high dimensional integrals. For a sufficient number of function samples, all the results 
agreed admirably with those obtained by routing simulations. 

Finally, we included a sleep scheme to demonstrate its effects on the local node density and the 
path dependence. For systems with a low p, the results imply that the independent model can be 
used, thus reducing computation time in calculating stochastic properties of the system. 
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.1 Derivation of Qi{ui+i) 

We outline how to calculate the rescaled mean measure for the region feasible region under the 
dependent model. The method is akin to calculating the equivalent area function Ai{ui+i) in the 
homogeneous case [H]. In fact, the derivation of Qi{ui+i) is included for completeness, and we 
refer the reader to previous work [M] for further details. 

We use the function A-!/)(u2) again to describe the angular width of the intersection of the source 
and current feasible regions. Subsequently, the rescaled mean measure on this intersection region 
is given by 



Qi\o(w2) = / A^{w2)dw2, 
which leads to the rescaled mean under the dependent model 

Ql(w2) = (5«i(m2) - (3l\o(u2) 




Figure 15: The form of intersection region depends on the U2 interval and the location of Xqi. 



It can be shown that the intersection angle expression is 

where Iq^ is an indicator function for when the intersection-point Xqi is below the baseline that 
runs from Xg to Xs (refer to Fig. [T5)) . and uqi is the sink distance of Xqi. On the second interval 
we obtain the intersection angle expression 

Aip{u2) = ^p^iu2) + i/'«i("2) - 6*1, "oi < U2 < "i- 
Recall under the independent model the rescaled mean integral 

r 

Q-jiu) = 2 / iljj{w)dw, 



where the angle function 



ip-yiu) — arccos 



, 

Thus, on the first interval, [i — r, uqi], we have the rescaled mean 



(3i\o('"2) =2 / 1p£{w2)dw2lx„^ 

Jl-r 

= Qt{u2)lx,n- 

On the second interval, [uqi,!*!], we have the slightly more complicated rescaled mean expression 

(■"2 

Ql\o('^2) = / ['4>l{w2)+ll^uAw2)-9i]dw2+Ql{uoi)lx^^^, 
Jl-r 

= ^ {Ql{u2) + Qui {U2) + 201 [moi - U2]) 
+ \ (Q£(uoi)[2Ii„, - 1] - OniKl)) • 
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This approach naturahy extends to the rescaled mean on the intersection of any two feasible regions. 
Consequently, for i > 1, under the one-hop dependent model the rescaled mean measure on the 
feasible region is given by 

Qi{Ux+l) = Qu,{Ul+l) - Qi\t-l{Ui+l). 



A Joint probability density 

We assume there is no sleep scheme, and note that to include a sleep scheme entails substituting 
the corresponding node density function p9p into the joint probability density. 

We consider the probability density of Qi conditioned on the event Ui = Ui. Under our inho- 
mogeneous Poisson model, the angle of any node is distributed uniformly around the sink (in the 
regions where nodes can exist). Hence, the conditional probability density under the dependent 
model 



yuA'^ll^l = Ul) = — — -, — :-, 



follows, and this expression also applies to the independent model as it has the same feasible region. 
We introduce the function ^^.{ui+i) to denote the total angular width of the feasible region given 
UiJ^i ~ Uj+i and the path Xi = Xi. The conditional probability density 

ITT \ ^V,iUt+l,di+l) 

5«,+Afi+iK^i+i — Ui+i) — — — — — , 

follows. Under the independent model the angular width function simplifies to 

The rescaled mean of the feasible region written as an integral 

J Ui — r 

gives the derivative of the rescaled mean measure 
The probability density 

follows. Hence, the joint probability density of the two random variables C/i+i and 6i+i is given by 

g,{ui+i,9i+i) = gs,iui+i)gu,+i{d'i.+i\Ui+i = u,+i,Xi = x{) 

where the spatially-dependent node density function has been introduced 

Ac, (wi+i, 6*^+1) = AIx),(ui+i,6ij+i). 
The joint probability density of the random variables Ui to [/„ and 0i to 9„ 

n 

9in-i){ui, ., u^, Oi, ., e„) = II \v,_, {u,,9,)e~^Q^-'(-"-\ (65) 

i=l 

follows, or in terms of hop advancements the equivalent expression 

n 



24 



